Distinct interacting cortical networks for stimulus-response and repetition-suppression

Non-invasive studies consider the initial neural stimulus response (SR) and repetition suppression (RS) – the decreased response to repeated sensory stimuli – as engaging the same neurons. That is, RS is a suppression of the SR. We challenge this conjecture using electrocorticographic (ECoG) recordings with high spatial resolution in ten patients listening to task-irrelevant trains of auditory stimuli. SR and RS were indexed by high-frequency activity (HFA) across temporal, parietal, and frontal cortices. HFASR and HFARS were temporally and spatially distinct, with HFARS emerging later than HFASR and showing only a limited spatial intersection with HFASR: most HFASR sites did not demonstrate HFARS, and HFARS was found where no HFASR could be recorded. β activity was enhanced in HFARS compared to HFASR cortical sites. θ activity was enhanced in HFASR compared to HFARS sites. Furthermore, HFASR sites propagated information to HFARS sites via transient θ:β phase-phase coupling. In contrast to predictive coding (PC) accounts our results indicate that HFASR and HFARS are functionally linked but have minimal spatial overlap. HFASR might enable stable and rapid perception of environmental stimuli across extended temporal intervals. In contrast HFARS might support efficient generation of an internal model based on stimulus history.

A ubiquitous finding in neuroscience is that neural responses to repeated stimuli are reduced compared to initial stimulus presentation, the phenomenon of repetition suppression (RS). RS has been shown in both single-unit studies in the monkey cortex 1-3 as well as noninvasive studies in humans using different techniques (for a review, see ref. 4 ). Several explanations for RS have been put forward, including adaptation or habituation, sharpening of representations, and reduction of prediction errors [4][5][6][7][8][9][10][11][12] . This reduction of responses to frequently occurring stimuli is associated with an enhanced response to unexpected events, establishing a mechanism for change detection 13,14 , with the probability of stimulus events accounting for a large proportion of neural variability 9 Most hypotheses on the mechanisms responsible for RS assume that what is suppressed is the stimulus-induced response. That is, the same neurons or networks that are initially responsive to the stimulus are the ones which are less active when the same stimulus repeats. Noninvasive studies in humans report that RS and stimulus response (SR) overlap, but these methods cannot distinguish nearby cortical activity.
A critical question remains whether RS is restricted to reducing the SR, in which case SR and RS should co-occur in the same electrodes (henceforward SR+RS+ sites), as suggested by scalp EEG and MEG recordings. Alternatively, SR and RS could be dissociated, but the circuits exhibiting SR and RS are intermingled and not resolvable by low-resolution scalp recording. This potential dissociation could be measured using a direct cortical recording of broadband high-frequency activity (HFA, 80-150 Hz), which is the key response frequency in previous ECoG (electrocorticography) studies [15][16][17][18][19] studying SR and RS.
Here, we utilized the high temporal and spatial resolution of direct cortical recordings from subdural ECoG electrodes to compare SR+ and RS+ signals in ten patients presented with trains of task-irrelevant auditory stimuli while attending a visual slide show to probe the automatic nature of RS. We show that while SR and RS both engage frontal, parietal, and temporal regions, they can be dissociated temporally and spatially in HFA. Critically, HFA SR+ and RS+ sites are distinctly modulated by θ and β low-frequency activity, respectively, with mutual information flow from SR+ to RS+ sites.

Results
For easy reference, the results sections correspond to the similarly enumerated sections of the methods section.
I -Stimulus response. We studied 412 channels (95 channels over frontal, 202 channels over temporal, and 115 channels over parietal regions across all subjects (see Fig. 1b)). Ninety-one channels showed a significant HFA modulation to auditory tones (of which 84 showed a stimulus response without RS, SR+RS-; Fig. 1b, c). Stimulus response occurred between 17-393 ms (SR max at 113 ms, p < 0.001) over multiple cortical regions (11 channels over frontal, 53 over temporal, and 20 over parietal regions, see Fig. 1b and Table 1), manifested as an increase in HFA power. None of the RS+ channels in our study showed a significant reduction of the HFA relative to the pre-stimulus period in response to the first standard.
II -Repetition suppression. RS (defined as both a significant F N standard and significant negative trend (r value) from S1 to S3) was found in 31 channels (of which 24 channels did not show SR, SR-RS+ Fig. 1d and Table 1) between 48-436 ms following stimulus onset (F max at 265 ms, p < 0.001, Fig. 1e). Channels designated as SR-RS+ showed no trend towards a stimulus response (average p = 0.41, SD = 0.1) even with a liberal, uncorrected significance criteria. We found three channels showing RE (both a significant F N standard and significant positive r value), in the parietal (N = 2) and the temporal (N = 1) cortex. None of these channels showed significant SR (i.e., they were SR-RE+ channels).
VI -SR+:RS+ integration. The distinct SR and RS cortical sites raise the question of whether SR+RS-and SR-RS+ sites interact. SR+RS-and SR-RS+ electrodes showed increased post-stimulus (116-270 ms) phase:phase coupling of θ SR+ and β RS+ (κ crit = 0.0043; κ max = 0.019; p < 0.00001, Fig. 3f) indicating significant interaction (Fig. 3g, h). Phase concentration coefficient κ differed significantly across frequency band pairs between 160-261 ms (F crit = 3.95 denotes the critical F value which the observed F values had to exceed to be considered significant; max F 2,288 = 8.27; p < 0.00001), due to stronger θ:β coupling (κ θ:β = 0.0091) than θ:α Fig. 1 Stimulus response and repetition suppression show different spatial profiles. a shows the auditory oddball paradigm. While the occurrence of deviants was unpredictable, S 1 , S 2 , and S 3 were always predictable. b shows the spatial distribution of HFA SR+RS-(green) and SR+RS+ (magenta) channels as marked dots against the background of a standard schematic brain using MNI coordinates that was also used for surgical planning. The remaining electrodes are marked by small white dots. c shows the HFA amplitude modulation of SR+ channels over time averaged across electrodes. The shaded area denotes the standard error across channels. d shows the spatial distribution of HFA SR-RS+ channels (blue) analog to b. Channels with magenta circles represent SR+RS+, as in b. e Modulation of repetition-suppression (F values) in SR-RS+ channels. The dashed blue line represents the significance threshold of F values. The shaded area denotes the standard error across the channels. coupling (κ θ:α = 0.0004; t = 3.25; p = 0.0016), and α:β coupling (κ α:β = 0.0009; t = 5.03; p < 0.0001 Fig. 3i, j).
VII -Information propagation. Between-sites mutual information (MI) analysis revealed early SR+ activity between 89-192 ms to be predictive of later RS+ activity between 190-226 ms (MI crit = 0.85 bits; MI max = 0.86 bits at 141 ms of SR + and 210 ms of RS+ time series; p < 0.00001. Fig. 3k) suggesting information propagation from SR+ to RS+ sites.

Discussion
Numerous studies report that the response to sensory stimuli decreases with repeating stimulation, a phenomenon known as repetition suppression (RS) or stimulus-specific adaptation 14 .
Noninvasive studies report substantial spatial overlap of stimulus response and repetition suppression, but such studies are limited in spatial resolution 20 . Thus, these methods are not well suited to examine whether RS reflects a stimulus response which gets reduced upon repeating stimulation, or might be a separate phenomenon of activity reduction relative to baseline in sites lacking stimulus response. This type of RS might reflect a shortterm memory mechanism, independent of stimulus response.
Here, we used intracranial EEG data in the context of repeating tones to measure the temporal, spatial, and spectral features of both phenomena. Unlike previous studies studying RS, we did not limit our analyses to channels showing stimulus response. If RS depended on sites being responsive to the stimuli, we should expect spatial overlap between the two, yet we found many sites showing exclusively SR and sites showing exclusively RS. In fact, a surprisingly small proportion of sites showed both SR and RS. Moreover, SR and RS exclusive sites also showed distinct spectral characteristics. SR sites showed higher θ power and θ:HFA phase-amplitude coupling (PAC) than sites with no SR, and RS sites showed stronger β power and stronger β:HFA PAC than sites showing no RS. Finally, we show that the two processes interact -the two types of sites show phase-phase coupling between their dominant theta and beta frequencies. Further, HFA SR peaks earlier than HFA RS effects, and HFA SR predicts HFA RS. Stimulus response was evident in all three lobes measured (frontal, parietal and temporal), manifested in high-frequency activity. All three regions also showed RS, with an amplitude decrease to repeated stimulation. Previous ECoG studies have No differences between channel sets were found within the α band (unsaturated color bars), gray dots represent single data. The right panel shows spectral power for the three consecutive standard trials separately, for the SR-RS+ channels (blue) and SR+RS-channels (green) separately. c, d HFA modulation by θ (c) and β (d) for SR-RS+ (blue) and SR+RS-(green) channels. e Kullback-Leibler divergence of θ and β for SR-RS+ and SR+RS-. f, g modulation of concentration index κ, indicating transient θ:β interaction, averaged across phases (f, the dashed horizontal red line shows the significance threshold) and for each θ phase (g). h θ-phase:β-phase coupling at κ max , gray lines represent pairwise combinations, red dashed line shows the average across all pairs. i the κ-comparison between θ:β, θ:α, and α:β (bottom: ANOVA F values). j Post hoc comparisons of θ:β, θ:α, and α:β, gray dots represent single data. (*** indicates p < 0.001) k Directed information propagation from SR+RS-to SR-RS+. In all panels, the shaded error margin around lines indicates the standard error across channels. shown adaptation of high-frequency activity, mostly located in the temporal and parietal cortex 21 , and here we show that lateral frontal cortex sites adapt to frequent auditory repetition as well, which has only been shown for low frequencies 22,23 .
What could be the mechanisms driving SR with no apparent RS? Under the traditional explanation of RS as a process of adaptation which recovers with time 24 , SR without RS might reflect sites with recovery times shorter than our inter-stimulus intervals, such that full recovery has been attained. Indeed, previous findings suggested variability in adaptation time constants 24 . The time constant of neuronal populations can be seen as reflecting the temporal resolution of the representation of the environment. Differences in resolution might allow distinguishing processing of coarse and fine-grained details 24 . Previous studies in vision show that the phase of θ activity across frontal and parietal regions is related to rhythmic sampling predicting visual detection performance 25,26 . Here, HFA SR+ modulated by θ activity (as evidenced by an increase in PAC), independent of RS, could be a mechanism to detect sensory evidence independent of context information or expectations.
An alternative explanation suggested by a reviewer is that SR-RS+ channels reflect the summed response of neurons excited by the stimulus and show repetition suppression (i.e. genuine SR+RS+) and of neurons that are a-priori inhibited by the stimulus. These inhibited neurons would "cancel out" the apparent SR in our mesoscopic recording of LFPs, making it look as if RS is present without SR. Without extensive single neuron measurements, it is difficult to rule out this possibility. However, we suggest that this possibility is less likely, considering the different spectral characteristics of SR+ sites and RS+ sites. If apparent SR-RS+ channels represent a composite signal of genuine SR +RS+ and inhibited neurons, then SR-RS+ channels should show comparable spectral signatures as other SR+ channels (namely SR+RS-and SR+RS+). However, both in the θ and the β range, SR-RS+ channels show a spectral signature which resembles other SR-channels (SR-RS-) and differs from SR+ channels. Specifically, in the θ range, SR-RS+ channels show θ power that is lower than in SR+RS-or SR+RS+ channels and is not different from channels neither responding to the stimulus nor showing repetition suppression (SR-RS-). In the β range, SR-RS+ channels also show different β activity compared to SR +RSchannels. Additionally, we found very few sites with a negative stimulus response to the first stimulus in a series, suggesting that inhibition is a rare response to the stimuli.
Within the classical notion of RS, it is harder to explain the finding of RS without an initial SR. One possibility is that rather than simple adaptation, repetition suppression without stimulus response reflects direct top-down inhibitory prediction signals. This is consistent with the conjecture that top-down prediction effects are carried by activity in the beta band. Previous studies argued that feedforward and feedback signals are distributed across cortical layers and segregated by spectral content 27,28 . They suggested that feedback signals target deep (infragranular) layers of the cortex with activity in the β band. In addition to enhanced β band activity, repetition suppression sites showed enhanced β:HFA PAC. Hence, repetition suppression embedded in β activity might reflect an internal model based on stimulus history. RS can be explained by neural sharpening 4,12,29,30 due to the fall-out of neurons that are not optimally tuned to stimulus features with repetition. Repeated sensory evidence strengthens intracortical inhibitory connections. This lateral interaction 31 may cause a decrease in the population response ('inhibitory sharpening´1 2 ). Even though we cannot directly test neural sharpening, lateral interactions can be an explanation for the SR +RSvs SR-RS+ sites. This sharpening may also be influenced by top-down inhibition as a component of hierarchical predictive coding.
Predictive coding (PC) schemes mostly assume that the stimulus response itself is suppressed when it is predicted, indicating that the evaluation of a stimulus likelihood precedes or occurs simultaneously with the bottom-up response to the stimulus, so that only deviations from prediction are registered. The current results show that RS effects followed SR in time, overlapping in time with only the late part of the stimulus response. This is consistent with prediction-error-like effects as the mismatch negativity (MMN), which overlaps in time the late part of the auditory N1 response from 100 ms after stimulus onset 32 . Thus, even if predictions are formed a-priori, they seem to affect mainly the later stages of processing.
The SR+ and RS+ ensembles showed distinct spectra and spatial and temporal layout, yet were not disconnected-activity in the θ and β range exhibited significant phase-phase coupling and HFA showed information transfer from SR+ to RS+ sites, which was local in time. The phase coupling signals originate from separable sites and were found exclusively between θ and β, suggesting that the θ:β cross-frequency phase coupling effects are not due to waveform shape 33 . Phase coupling might enable temporally precise coordination of neuronal processing by establishing systematic spike-timing relationships among functionally distinct oscillatory assemblies enabling functional integration and coordination 34,35 . The functional meaning of this coupling in the present case remains to be investigated.
Another sign of communication between the sites is that HFA in the SR+ sites reduced the uncertainty of HFA responses in the RS+ sites. MI is based on information theory 36 and formally characterizes the information content of neural responses and interactions between these responses. In previous studies, MI has been applied to multiunit recordings and local field potentials in nonhuman primates 25,[37][38][39] and intracranial data in human 40 . Importantly, MI makes no assumptions about the content of the signal itself but only that it changes as a function of time. The strength of this approach is that it characterizes the nonlinear relationship between two different neural responses. We found a clear exchange of information from SR+ to RS+ sites but not vice versa. Under the assumption that the repetition suppression reflects the process of top-down model predictions, this result supports the notion that stimulus responses (or prediction errors under the PC framework 41 ) inform the generation of the internal model. Taking together the finding of phase:phase coupling between θ and β activity, the fact that each band respectively Table 2 Comparison of PSDs between channel types SR+RS-, SR-RS+, SR-RS-, and SR+RS+ channels, respectively, each for the three frequency bands θ, β, and α based on the interaction between both factors. modulated the local HFA amplitude in its region (PAC), and the finding of directional MI between of SR+ to RS+ sites, supports the proposal of cross-talk between bottom-up (θ) and top-down (β) effects.
Limitations. Repetition effects were investigated in depth in previous studies in vision and audition across different species using a large repertoire of recording techniques 19 . These studies are often motivated by behavioral repetition priming, showing that repetition leads to improved identification of stimuli 42 . Repetition-suppression is assumed to reflect statistical learning 43,44 and contributes to sensory memory update 45 by tracking stimulus history 46 . Most repetition suppression studies compared responses between a first and second presentation which show the strongest repetition suppression effects in noninvasive recordings 20,[47][48][49] . Here, we calculated RS across the first three standards since we conjectured a monotonical decrease 50 until the fourth standard 45,47 . Several previous studies showed local PAC as in our study (phase of low frequency and amplitude of higher frequency of the same broadband signal) across species and different recording techniques [51][52][53][54][55][56][57][58][59][60][61][62][63][64] . However, since PAC supports information processing, genuine coupling in contrast to spurious correlations must be shown. Spurious coupling can result from filtering non-sinusoidal signal 34 creating artificial coupling at distinct frequencies, especially in local PAC. However, commensurate with phase:phase coupling, PAC is more likely genuine if the higher frequency is exclusively coupled to only one of two distinct low frequencies originating from separable neuronal processes. We indeed found a double dissociation of HFA SR+ coupled to θ but not β and vice versa for HFA RS+ in the temporal interval of coupling between θ and β. This differential coupling finding provides evidence for two distinct processes 34 .
In our study, we focused on the repeated and hence predictable part of the auditory stimulation. How repetition and expectation suppression interact is still debated. Our paradigm does not allow us to conclusively dissociate between repetition independent of expectation. However, recent evidence 65 suggests that subjects are less likely to apply overt expectations (even when available) when the stimuli are not task-related (as was the case in our study), suggesting that when attention is directed elsewhere, expectation vanished despite robust repetition suppression being still evident 65,66 .

Conclusion
Our critical finding is that SR and RS dynamics were temporally and spatially distinct. Our results highlight the role of distinct processes in computing a stimulus response and feedback signals which are functionally linked but do not completely overlap.

Methods
Patients. Ten patients (mean age 32, SD = 9.74) undergoing pre-surgical monitoring for drug-resistant epilepsy 67 with subdural electrodes participated in the experiment after providing their written informed consent. Recordings took place at the University of California, San Francisco (UCSF) (five patients) and at the Dept. of Epileptology (Krankenhaus Mara), Bielefeld University (five patients) and were approved by the local ethics committees. Data from these patients were preprocessed in an analogous manner as reported 18 .
Stimuli. Participants listened to stimuli consisting of 180 ms long (10 ms rise and fall time) harmonic sounds with a fundamental frequency of 500 or 550 Hz and the 3 first harmonics with descending amplitudes (−6, −9, −12 dB relative to the fundamental). The stimuli were generated using Cool Edit 2000 software (Syntrillium, USA). They were presented from loudspeakers positioned at the foot of the subject's bed at a comfortable loudness.
Procedure. While reclined in their hospital bed, participants watched an engaging slide show while sound trains were played in the background. Sound trains included high probability standards (p = 0.8; f 0 = 500 Hz) mixed with low probability deviants (p = 0.2; f 0 = 550 Hz) in blocks of 400 sounds, with a stimulus onset asynchrony (SOA) of 600 ms. Hence, in each block, 320 standard tones and 80 deviant tones were presented. In different blocks, the order of the sounds was either pseudorandom, with a minimum of three standard tones before a deviant (irregular condition), or regular, such that the standard stimulus was repeated exactly four times before a deviant was presented (S-S-S-S-D-S-S-S-S-D-…, Fig. 1a). Thus, in the regular condition, the fourth standard tone was fully predictable, whereas in the irregular condition, the fourth stimulus could be either a standard or a deviant, and prediction was not possible. The current report examines the responses to the repeating standards. The deviance-related responses from the subset of the subjects recorded at UCSF was previously reported 18 .
Data recording and preprocessing. ECoG was recorded at UCSF using electrode grids equipped with 64 platinum-iridium-electrodes, arranged in an 8 × 8 array with 10 mm center-to-center spacing (Ad-Tech Medical Instrument Corporation, Racine, Wisconsin). At The Mara, Bielefeld, ECoG was recorded via electrode strips (single strips or parallel arrangement of strips; white dots in Fig. 1b, d represent all electrode locations) using a Nihon Kohden amplifier (Tokyo, Japan). Electrodes were positioned based solely on clinical needs. The exposed electrode diameter was 2.3 mm. The data at UCSF were recorded continuously throughout the task at a sampling rate of 2003 Hz. At The Mara, the sampling rate was 2000 Hz in the case of four subjects and 1000 Hz in one subject. We used Matlab 2013b (Mathworks, Natick, USA) for all offline data processing. After visual inspection, we excluded channels exhibiting ictal activity or excessive noise from further analysis. In the remaining "good" channels (see Table 1), we then excluded time intervals containing artifactual signal distortions such as signal steps or pulses by visual inspection. Finally, we re-referenced the remaining electrode time series by subtracting the common average reference x c ðtÞ ð1Þ calculated over the n good channels from each channel time series x c . The resulting time series were used to characterize brain dynamics of responses to repeated auditory stimulus presentation. For high-frequency signals, we band-pass filtered each electrode's time series in the high-frequency range (80-150 Hz). All filtering was done with zero-shift infinite impulse response (IIR) filters [Butterworth filter of order 4: filtfilt() function in matlab]. We obtained HFA by calculating the analytic amplitude A f ðtÞ by Hilbert-transforming the filtered time series. We smoothed the HFA amplitude time series such that the amplitude value at each time point t is the mean of 10 ms around each time point t. Filtering was done for each trial (-1 s to 2 s around stimulus onset-sufficiently long to prevent any edge effects during filtering).
Data analysis. We conducted the following analysis steps explained in detail below. First, we defined stimulus response of the HFA (I-Stimulus-responsive activity modulation). Next, we parameterized response attenuation of cortical HFA responses to repeated standard tones using a time-resolved ANOVA (II -Repetition-suppression). We then compared the temporal and spatial profile of stimulus response and repetition suppression (III -Comparison of stimulus response and repetition suppression). Next, we tested to what degree the HFA SR and RS are associated with low-frequency activity (IV -Comparison of low-frequency specific modulation). We then tested whether HFA SR and RS were modulated by distinct neural populations in low frequencies (V -Cross-frequency modulation) and for functional integration between low-frequency neural populations (VI -SR:RS integration). Finally, we assessed information flow between SR + and RS + channels using time-resolved mutual information to test for the directionality of information propagation (VII -Information propagation).
I -Stimulus-responsive activity modulation. We identified stimulus-responsive channels SR+ showing a significant HFA modulation following the onset of standard stimuli using the following steps. We first averaged stimulus-locked HFA responses across all standard trials. To apply a sensitive measure which takes into account that HFA responses can occur with a delay and/or can be transient, we calculated A HFA as averaged activity modulation across five different intervals, each with a duration of 100 ms following the stimulus onset (starting at 0, 50, 100, 150, and 200 ms). This allowed us to select fast and delayed responding channels. We then calculated the average baseline activity B HFA across the 100 ms preceding the stimulus onset. For the stimulus-related HFA, we subtracted baseline B HFA from the activity modulation A HFA following stimulus onset for each of the five temporal intervals. The difference between B HFA and A HFA was compared against a surrogate distribution derived from randomly shifted time series (1000 permutations). In each iteration, the time series of each channel and each trial were shifted (circular shift of the entire trial time series between −0.1 and 0.6 s) separately, and new (surrogate) trial averages ( B surr and A surr ) were calculated from the shifted trials, resulting in a surrogate distribution of differences between baseline and activity. The comparison of observed difference values with the surrogate distribution results in a p value for each channel and time interval. To control for multiple comparisons, we corrected the p values by applying the false discovery rate (FDR) 68 method across all channels and time intervals. Channels with a q < 0.05 (where q is the false discovery rate) in any of the five intervals were classified as showing a significant HFA modulation following the standard stimuli and were denoted as SR+, whereas the remaining channels were labeled as SR-. To determine the amount of evidence for a change over baseline, we compared HFA values at each time point with HFA values in the baseline interval separately for each channel across trials, using Bayes factor (BF; bf.m toolbox in MATLAB http://klabhub. github.io/bayesFactor/). BF >3 is considered strong evidence for a difference (the difference is three times more likely than no difference) and BF <1/3 supports null effects 69,70 .
II -Repetition-suppression. RS is defined as attenuated amplitude to repeated stimulus presentation. Hence, this definition is twofold: (i) a change of amplitude which is (ii) monotonically decreasing with the number of repeated stimulus presentations. While (i) refers to statistical differences in brain response with repetitions, (ii) assumes a specific model of response attenuation. We thus grouped trials according to the number of standards in a train in three groups (S 1 , S 2 , and S 3 ) since only the first three standards in a train can be expected in both conditions. To parameterize the amplitude modulation with stimulus repetition (i), we ran a one-way ANOVA with a factor number of standards for each electrode (with trials as a random variable), regardless of whether it was SR+ or SR−, at every time point, both in the regular and irregular condition. This leads to an F value time series (main effect: F N standard ) for each channel in each condition. Significant F values only define differences between numbers of preceding standards but not the exact model of monotonical decrease of neural responses. We tested (ii) the model of a monotonic neural amplitude decrease across the number of repetitions. For each electrode, at each time point, we calculated the Pearson correlation between the mean HFA and the number of sequential preceding standards, yielding a time series of correlation coefficients (r) for each channel. We compared each F value and r value against a surrogate distribution constructed under the null assumption of no difference or no correlation, respectively. This surrogate distribution was constructed by randomly reassigning the labels (S 1 , S 2 , and S 3 ) to the single trials in 1000 permutations for each channel. This leads to 1000 surrogate F Nstandard and r value time series. We assigned a p value to each F and r value within the surrogate distributions. The p values for F and r were corrected for multiple testing by applying the FDR. F and r with q < 0.05 were classified as significant. Finally, channels showing intervals with the conjunction of both significant F N standard and significant negative r value were considered as showing significant RS and were labeled RS+. In contrast, channels with temporal intervals of significant F value and a significant positive r value show repetition enhancement (RE+). In addition, to ensure that RS+ shows an S1 > S2 > S3 relationship, we computed the amplitude decrease over a train of tones across these channels. We averaged mean responses following S 1 , S 2 , and S 3 across RS+ channels. We then calculated the differences between those responses (response following S3 vs response following S 2 and response following S 2 vs response following S 1 , respectively) at each time point and summed the sign function of differences (−1 and +1 for negative and positive differences, respectively).
III -Comparison of stimulus response and repetition suppression. We tested the spatial and temporal relation between stimulus response (channels selected in step I) and repetition suppression (channels selected in step II) in the following way. If SR and RS are multiplicatively related 21 , channels with stronger stimulus response should show stronger repetition suppression. Alternatively, suppression could be subtractive (the same amplitude reduction regardless of SR), in which case RS will be fixed, that is, not dependent on SR. To examine this, we tested the correlation between SR and RS measures across channels. For each channel, we averaged HFA in response to S 1 in the time window of strongest SR, averaged F N standard ( F) and r ( r) values separately (section II) in the interval of significant RS, and separately correlated the S 1 HFA with F and with r. This correlation was tested both across all channels and across channels showing a significant stimulus response (as defined in I above) and/or a significant RS+ (as defined in II above).
IV -Comparison of dominant band power. We estimated the power spectral density (PSD) in each trial (collapsing across all S 1 , S 2 , and S 3 trials) separately for SR+ and RS+ channels using Welch's method based on the FFT 71 . Specifically, for each channel, we calculated PSD as a function of frequency (1-30 Hz, 1 Hz steps) in each trial in temporal intervals of 100 ms (due to the short SOA of 600 ms) between −0.1 s to 0.3 s in steps of 50 ms. The resulting PSD values were averaged across trials yielding one PSD for each channel. We then averaged across three canonical low-frequency bands θ (4-8 Hz), α (8-12 Hz), and β (15-30 Hz) and compared the resulting power estimates across channels using a two-way ANOVA with the factors channel type (SR+RS-, SR-RS+, SR-RS-, and SR+RS+) and frequency band (θ, α, β). Post hoc we compared the power estimates between channel sets, separately for the three different frequency bands by computing t values and comparing those to the critical t crit . t crit denotes the critical t value which the observed t values had to exceed to be considered significant. To correct for multiple comparisons, p values were assigned to each t value within a surrogate distribution constructed by randomly assigning labels (SR+RS-, SR-RS+, SR-RS-, and SR+RS+ individual channels) and corrected by applying the FDR. The adjusted p values are labeled q. T values with a corresponding q < 0.05 (corrected p value) were classified as statistically significant.
V -Cross-frequency coupling. The interplay between activity at distinct frequencies is proposed to be regulated via cross-frequency phase coupling (CFPC 72 ) and via phase-amplitude cross-frequency coupling (PAC; see below [73][74][75] . Phase-amplitude cross-frequency coupling (PAC) is a mechanism that has been proposed to coordinate the timing of neuronal firing within local neural networks (see ref. 73 for a review). We utilized conventional cross-frequency coupling metrics 51,76 to test for differences in coupling of HFA to low-frequency bands in SR+ vs. RS+ channels. We calculated the instantaneous phase for low-frequency activity (see below) for each SR+ and RS+ channel time series. In the temporal interval of coupling between low-frequency networks as determined in the previous step IV, we divided both the θ and β cycle separately in 50 equally spaced bins ranging from -π to π and computed the average HFA for all trials within a 45-degree window centered on every phase bin 25 . The resulting HFA histograms-each containing 50 valueswere averaged separately for SR+ and RS+, separately for θ and β activity. We then calculated the normalized Kullback-Leibler divergence (KLD) of the observed distribution against a uniform distribution to quantify how strongly the observed distribution of HFA SR+ and HFA RS+ were modulated by the phase of the θ or β activity. The obtained KLD were compared in a two-way ANOVA with the factor channel type (RS+ vs. SR) and frequency band (θ vs. β). The interaction effect of the ANOVA describes the double dissociation of HFA of SR+ and RS+ PAC to θ and β networks, respectively.
VI -SR:RS integration. We hypothesized that while the amplitude of HFA is modulated by the phase of the low frequencies within a population (PAC 14 ), communication across populations in low frequencies might be achieved via CFPC 34,52,77 . CFPC is defined by a nonrandom phase difference between oscillations, enabling temporally precise coordination or integration among functionally distinct oscillatory networks 34 . We tested whether phases of canonical lowfrequency bands f 1 and f 2 -either θ:α, α:β, or θ:β-are aligned and whether this interaction is modulated in time. To that end, we calculated the instantaneous phase of SR+ and RS+ channel time series and binned phase time series in intervals of one cycle of the f 1 (133 ms for θ and 100 ms for α) centered on time points between −200 to 300 ms following stimulus onset in each trial. In each temporal interval in each trial, we divided the f 1 cycle into 50 equally spaced bins ranging from -π to π and registered the f 2 phase at each bin. This was done for each pair of SR+ and RS+ channels, excluding SR+RS+ channels within each single recording session and separately for each patient, in which we found both RS and SR channels. This results in an f 2 phase angle distribution at each f 1 phase and at each time point for each pair. For each distribution, we calculated the concentration coefficient κ (reciprocal value to variance) across all f 2 phase angles at each f 1 phase at each time point. κ time series of each SR+/RS+ channel pair were baseline corrected by subtracting the mean κ values in the 200 ms preceding stimulus onset. We then averaged κ values across all f1 phases leading to a κ time series for each pair of SR+/RS+ channels. Temporal intervals of high κ (low variance of f 2 phases) indicate coupling of the f 2 phase to f 1 phases. Low κ, on the other hand, indicates that f 2 and f 1 phases are unrelated. Each κ-value was compared against a surrogate distribution. In 1000 runs, we shifted phase time series in each trial and each channel separately and calculated surrogate κ values.
To correct for multiple comparisons, p values were assigned to each κ-value within the surrogate distribution and corrected by applying the FDR procedure. We then compared κ θ:β with κ θ:α and κ α:β time series. To parameterize the difference in f 1 :f 2 coupling, we ran a one-way ANOVA with factor frequency pairs (θ:β, θ:α, and α:β) at each time point. This leads to an F value time series representing the difference in coupling strength between θ, α, and β pairs. Each F value was compared against a surrogate distribution and a p value within the surrogate distribution was assigned. This surrogate distribution was constructed by randomly reassigning the labels (θ:β, θ:α, and α:β) to the f1:f2 time series in 1000 permutations leading to 1000 surrogate F values. To correct for multiple comparisons, we corrected the p values resulting from comparing F values against the surrogate distribution by applying the FDR.
VII -Information propagation. Using mutual information (MI), we tested whether there is directed information flow between SR+ and RS+ channels. MI measures how much a random variable can be predicted by another random variable. Specifically, MI quantifies the uncertainty about one random variable given knowledge of another variable and is given in units of bits. Mathematically MI was calculated as the sum of entropies of SR+ and RS+ minus their joint entropy. In each subject, we calculated MI between all pairs of RS+ and SR+ channels, excluding SR+RS+ channels within each subject. This was done Where H(SR) and H(RS) stands for the entropy of the SR+RS-and the SR-RS + channel, respectively, and H(RS, SR) designates their joint entropy. The denominator standardizes each value by the maximal achievable information value. We iterated through all intervals around each time point of RS+ and SR+ channels resulting in a matrix of MI values quantifying which temporal interval of the HFA time series of SR+ channels predicts the HFA time series of RS+ channels and vice versa for each pair of channels. We then averaged the MI across all pairs of channels. We then compared each MI-value against a surrogate distribution. This surrogate distribution was constructed by randomly shifting SR+ and RS+ time series of single channels and averaging across subjects in 1000 permutations. In each run, we repeated the same analysis as outlined above. This leads to 1000 surrogate MI values. The resulting p values for the MI values relative to the surrogate distribution were corrected by applying the FDR.
Reporting summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Data availability
The datasets generated and/or analyzed during the current study are available in the Open Science Foundation repository (https://osf.io/ceutw/?view_only=1d6ba767ebe3458 ca26ade4945972d8c).